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THEORY  AND  APPLICATION  OF  THE  WAVELET 
TRANSFORM  TO  SIGNAL  PROCESSING 


1.  INTRODUCTION 

Most  scientists  and  engineers  are  familar  with  the  Fourier  representation  of  a  func¬ 
tion.  In  most  cases  of  practical  interest,  it  is  well  known  that  a  piecewise  continuous 
periodic  function  gr{t)  can  be  equated  almost  everywhere  with  a  weighted  countably  in¬ 
finite  sum  of  exponential  functions,  i.e., 

Sr(t)=-  t  (1) 

n=~oo 

where  gn  is  the  nth  Fourier  coefficient,  and  is  given  by 

=  (2) 

This  is  the  classic  Fourier  series.  The  Fourier  representation  for  a  piecewise  continuous 
nonperiodic  function  g{t)  has  the  integral  formulation 

g(t)  ='  /”  GUVlw,,df,  (3) 

J  —  OO 

where  G(f )  is  the  Fourier  transform  of  5(f),  and  is  given  by 

<?(/)  =  r  (4) 

J  —  OO 

The  development  of  Fourier  theory  constitutes  a  significant  portion  of  classical  mathemat¬ 
ical  analysis.  It  has  also  been  a  valuable  tool  in  physics  and  engineering  by  contributing  to 
the  solution  of  problems  in  linear  system  theory,  thermodynamics,  and  quar'  urn  physics, 
to  name  a  few. 

Despite  their  wide  use,  some  problems  arise  in  the  interpretation  of  the  classic  Fourier 
representations.  For  example,  if  g(t)  is  a  function  that  is  nonzero  only  in  an  interval  (it 
has  compact  support),  then  the  Fourier  transform  implies  that  this  time-limited  function 
is  a  summation  of  complex  exponential  functions,  each  having  support  over  the  entire 
real  line.  Furthermore,  one  cannot  associate  features  of  the  time  function  g(t)  with  any 
specific  value,  or  range  of  values,  of  /  of  the  transform  G(f).  In  other  words,  the  transform 
exhibits  no  locality  of  time;  a  transient  feature  in  g(t)  contributes  to  G(f)  at  all  values  of 
/• 
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In  light  of  these  observations,  researchers  have  sought  alternatives  to  the  Fourier 
representations  that  display  some  sense  of  locality  in  time.  It  was  recently  discovered  that 
certain  functions,  those  whose  amplitude  is  significant  over  a  finite  interval,  can  serve  as  a 
transform  kernel  for  square  summable  functions,  provided  they  obey  a  regularity  condition 
[1-5].  These  functions  are  defined  as  a  time-dilated  and  time-shifted  version  of  a  specified 
function  i>(t)  called  a  wavelet.  Thus,  the  transform  kernel  is  ip( s (t  —  r)),  where  s  E  [0,  oo), 
and  r  6  (  —  00,00).  If  we  integrate  the  function  g(t)  with  respect  to  this  kernel,  the  result 
is  the  wavelet  transform  4>b(s,t),  which  is  analogous  to  the  Fourier  transform  since  it  has 
an  integral  representation  for  both  the  forward  and  inverse  transform.  There  is  also  a 
version  of  it  analogous  to  the  Fourier  series  in  that  it  represents  a  function  as  a  countable 
sum  of  wavelets. 

Many  uses  for  the  wavelet  transform  have  been  proposed  [1].  Among  them  are 
the  identification  and  localization  of  transients  in  time  series  and  the  compact  coding 
of  images  if  the  concept  of  a  wavelet  is  extended  to  two  dimensions.  In  this  report  we 
are  concerned  with  the  application  of  the  wavelet  transform  in  signal  theory  and  signal 
processing.  The  first  part  of  the  report  presents  several  theorems,  some  of  which  are 
results  that  have  either  been  left  unproven  in  other  references  or  have  only  been  alluded 
to  by  other  authors.  Others  are  new,  and  present  a  deliberate  attempt  to  reformulate 
some  standard  signal  theory  in  the  context  of  the  wavelet  transform.  The  last  part  of  the 
report  presents  the  application  of  the  wavelet  transform  to  two  common  signal  processing 
problems:  filtering  and  deconvolution. 

Throughout  this  report  we  maintain  the  following  notational  conventions:  time  sig¬ 
nals  are  represented  by  lower  case  italic  as  in  g{t),  and  their  associated  Fourier  transforms 
are  always  written  with  the  same  letter  but  in  upper  case  italic  as  in  G(f).  Vectors  are 
written  in  bold  lower  case  italic  as  in  as,  and  matrices  are  written  in  bold  upper  case  italic 
as  in  A.  Finally,  it  is  assumed  that  the  reader  has  some  familiarity  with  real  analysis. 


2.  THE  WAVELET  TRANSFORM 


As  the  Fourier  transform  does,  the  wavelet  transform  also  provides  a  representation 
of  the  elements  in  the  set  of  all  square  summable  functions.  This  is  the  space  L2(R ), 
where  R  denotes  the  real  line.  The  wavelet  transform  of  the  function  g(t)  is  given  by 

w  {5(0}  =  MS1T)  =  V*  f_  T))dt,  (5) 

where  ip(t)  is  referred  to  as  the  analyzing  wavelet,  s  is  the  dilation  variable,  and  r  is  the 
delay  variable.  This  defines  a  mapping  from  the  one- dimensional  (1-D)  time  domain,  to 
the  two-dimensional  (2-D)  dilation/delay  space  defined  by  the  (s,r)  plane.  The  factor  y/s 
appears  because 

3  I  \if>{s(t  -  r))\2dt  =  [  \i/>{t)\2dt,  (6) 

j  —  00  J  —  OO 
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thereby  acting  as  a  normalization  factor.  This  implies  the  relative  amplitudes,  or  powers 
(square  magnitude),  of  the  wavelet  transform  at  two  different  points  in  the  (s,t)  plane 
can  be  compared  meaningfully. 

Note  in  the  mathematical  literature  that  the  time-dilated  and  time-shifted  version 
of  the  analyzing  wavelet  is  defined  as 


where  a  —  1  /s,  and  a  €  [0,oo).  We  choose  the  definition  by  using  s  rather  than  a  since 
it  is  often  interpreted  as  the  time-dilation  due  to  the  Doppler  effect  caused  by  a  moving 
point  reflector.  That  is,  if  a  transmitter  is  stationary  and  a  reflector  is  moving  with  a 
positive  radial  velocity  of  v ,  then 

1  —  vie  2v  .  . 

3  =  - - —  ~  1 - ,  (8) 

1+v/c  c 

where  c  is  the  propagation  speed.  Thus,  if  we  transmit  the  signal  g(t),  then  we  receive 
y/sg(st).  Adopting  the  convention  of  modeling  time-dilation  by  the  variable  s  will  make 
<the  results  presented  in  this  report  immediately  applicable  to  those  working  in  the  field 
of  wideband  echo-location  systems. 


Given  the  definition  of  a  forward  transform  in  Eq.  (5),  we  naturally  seek  a  formula 
for  transforming  from  the  (s,t)  plane  back  to  the  time  domain.  If  the  wavelet’s  Fourier 
transform  '?(/)  meets  a  regularity  condition,  namely 


mn  i5 


df  <  oo, 


and  |$(/)|  =  |$(-/)|,  then  an  inversion  formula  for  the  wavelet  transform  exists,  and 
g(t)  can  be  recovered  from  <^fl(s,T)  through  the  formula 

5(0  =  77-/  /  \/s  <t>g(3,r)il>(s(t- T))drds.  (10) 

J * =0  J r=—oo 

For  completeness,  a  proof  of  this  result  is  given  below  since  most  published  versions  of  it 
can  only  be  found  in  obscure  references. 

Proof  of  the  Inversion  Formula  -  We  begin  by  substituting  Eq.  (5)  into  Eq.  (10)  to 
get 

1  =  ir  /  /  s  f  g{y)ip*{s{y -r))dy  ip{s(t-T))dTds.  (11) 

J  a— 0  J T—  —  oo  «/  —  oo 

However,  it  can  be  shown  that  the  Fourier  transform  of  ip(s(t  —  r))  is  given  by 

T  |tf(*(<  -  r))|  =  J°°  -  T))e~fl*}tdt  =  -  (12) 
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Therefore,  by  invoking  Parseval’s  Theorem,  using  Eq.  (12),  and  defining  G(f )  to  be  the 
Fourier  transform  of  g(t )  ( G(f )  exits  since  g(t)  £  L2(R)),  one  finds  that 

i  =  —  r  r  s  r  <?(/)**■— K^df  r  drds 

Cfj,  J f =0  Jt=  —  oo  J—oo  3  J—oo  3 


=  t ~r  r  s-iG{f)v'{fis)\r  r *(U3y2^+«t-T»dTd( d/ds.  (13) 

J §=0  J  }  =  - oo  |/-oo  J  —  oo 

The  bracketed  factor  in  the  integrand  of  the  last  line  of  Eq.  (13)  can  be  reduced  as  follows: 

/“  f°°  *(£/a)ea2,r</T+*(f~T))elT#  =  /°°  [  r  e2jKiTdr 

J —OO  J —  OO  J  —  OO  I «/  — oo 

=  f  sx()(s(t  —  T))e32yr^rdr 

J  —oo 

=  r  ^{zy^l^dz  ej2ir/t 

J —oo 

=  tf(//s)ej2,r/t.  (14) 

Thus,  substituting  Eq.  (14)  into  Eq.  (13)  yields 

1  =  1  T  r  s-G(/)|»(//i)|Je"",rf/^ 

J  t=0  J  f=—oo 

=  <15) 


Concerning  the  bracketed  factor  in  the  integrand  of  Eq.  (15),  if  /  /  0,  then  the  change  of 
variable  77  =  3/  f  yields 

r  ,-'\<i(f/s)\2d,  =  r  =  (is) 

Jo  Jo  TJ 

Note  that  this  is  true  regardless  of  the  sign  of  /  since  |\P(/)|  =  |®(— /)|  by  hypothesis. 
We  now  consider  the  value  of  the  bracketed  factor  as  /  — >  0.  In  this  case 

lim  f  /s)\2ds  =  lim  lim  /  s_1|\P(//a)|2da 

/  — *0  Jo  /-> 0  o->0+  Ja 

///° 

=  lim  lim  /  z  1|'I'(z)|2d.z 
/— *0[a— .0+  Jfa 

=  lim 
0 

=  Ct  (17) 

This  result  can  also  be  deduced  by  noting  that  for  any  |/|  >  0,  however  small,  Eq.  (16) 
always  holds.  Hence,  the  value  of  the  limit  of  the  integral  in  Eq.  (16)  is  C Regardless 
of  its  value  at  /  =  0,  we  at  least  know  that 


f  \^f(f /s)\2ds  —  for  almost  all  /. 
Jo 
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In  light  of  Eq.  (18),  Eq.  (15)  reduces  to 

/=  l°°  G(f)e*'»df^g(t),  (19) 

J  ~oo 

which  proves  the  inversion  formula  in  Eq.  (10).  This  completes  the  proof. 


3.  PROPERTIES  OF  THE  WAVELET  TRANSFORM 

We  now  state  and  prove  three  properties  of  the  wavelet  transform.  The  first  property 
concerns  the  wavelet  transform  of  a  time-dilated  and  time-shifted  signal. 

Theorem  1:  Let  (f> fl(s,r)  be  the  wavelet  transform  of  the  function  g(t).  Then  for 
•So  >  0, 

W  {v/^(30(<  -  r°))}  =  ~  T°))‘  (20) 

Proof  of  Theorem  1:  Equation  (20)  follows  directly  from  the  definition  of  the  wavelet 
transform  by  substituting  y/s^g(s0(t  —  r))  for  g(t)  in  the  integrand  of  Eq.  (5),  and  by 
invoking  the  change  of  variable  y  =  s0(t  —  r0).  This  completes  the  proof. 

The  theorem  above  shows  how  distortion  in  the  time  domain  affects  the  wavelet 
transform.  For  example,  if  s0  >  1,  we  see  that  the  support  of  the  wavelet  transform  in  the 
(a,r)  becomes  wider  in  s  and  narrower  in  r.  Delaying  a  function  in  time  merely  delays  its 
wavelet  transform  along  the  r  axis. 

As  in  the  case  of  the  Fourier  transform,  one  is  often  concerned  about  how  rapidly 
the  transform  G(f)  decays  as  |/|  — >  0.  This  ‘decay  rate’  is  related  to  the  continuity  of 
g(t).  The  following  theorem  describes  the  decay  rate  of  the  wavelet  transform  (j>g(s,r)  at 
points  of  continuity,  and  at  points  of  a  jump  discontinuity  for  the  function  g{t ),  providing 
g(t )  is  of  bounded  variation. 

Theorem  2:  Let  g(t)  6  L2(R)  and  ^( t )  €  Ll(R)  f|  L2(R).  Suppose  g(t)  is  of  bounded 
variation  in  a  neighborhood  of  to,  then  the  wavelet  transform  of  g(t)  with  respect  to  the 
wavelet  has  the  property 

(i)  <f>g(s,t0)  — ♦  0  as  s  — >  oo. 

Furthermore,  if  for  all  8  >  0, 

s  I  \x(;(z)\2dz  — *  0  as  s  — >  oo,  (21) 

then  if  g(t)  has  a  finite  jump  discontinuity  at  t0 ,  then 

(*'*)  \/^<Pg{si  ^o)  *  g{t o  +  g{to  )V*+  as  3  — ►  oo, 
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where 

V>+  =  / 

Jo 

=  f  xj}(t)dt.  (22) 

If  <7(£)  is  continuous  at  to,  then 

(Hi)  y/s<f>g(s,t0 )  — >  0  as  s  — >  oo. 


Proof  of  Theorem  2:  We  note  that  for  any  8  >  0, 

4>g(s,t 0)  =  /  aiWisit  -  to))dt 

J  —  oo 

=  [  +  /  y/sg(t)ip*(s(t  -  t0))dt 

--  Il  +  1 2  ■ 

Consider  /i  and  /2,  one  at  a  time. 


(23) 


We  first  show  that  A  — >  0  as  s  — »  oo.  By  the  Schwartz  inequality 

«<IAI!  <  Ml/  <W<(»  -  <o))l!* 

=  W/,  Mz)\ 2dz.  (24) 

But  since  8  >  0, 

lim  f  \if)(z)\2dz  =  lim  f  +  f  \if(z)\2dz  =  /  +  [  |^>(z)|2dz  =  0  +  0  =  0.  (25) 

1 >00«/|2|>#$  •—>ooJ_  oo  J  sS  %/— OO  Jo O 

(This  could  also  have  been  proved  by  using  the  density  of  the  step  functions  in  L2(R).) 
Thus,  by  applying  the  squeeze  theorem  for  limits  to  Eq.  (24)  one  finds  that 

lim,  |A|2  =  0  =>  lim  I\  -  0.  (26) 

s — *oo  a—*  oo 

Furthermore,  if  Eq.  (21)  is  true,  then  by  a  similar  development 

0  <  |\/*A|2  <  ||<7||2  /  s\ip(z)\2dz  — »  0  as  s  — *  oo,  (27) 

which  proves 

y/sli  — *  0  as  s  — ♦  oo.  (28) 

Now  consider  A-  It  is  sufficient  to  consider  the  case  where  g(t)  and  ift(t)  are  real 
valued.  Because  g(t)  is  of  bounded  variation  in  a  neighborhood  of  t0,  say  [£0  —  8,  to  +  5],  we 
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know  g(t)  is  the  difference  of  two  increasing  function  (that  are  also  of  bounded  variation). 
Therefore,  it  is  also  sufficient  to  consider  the  case  of  g{t )  increasing,  and  with  a  finite 
jump  discontinuity  at  t0. 


Using  the  second  mean  value  theorem  for  integrals,  we  note  that  for  some 
'•?  €  (toi  to  +  £]) 


/•to+i  rrj  fto+o 

/  g(t)sip(s(t  -  t0))dt  =  +  g(t)sip(s(t  -  t0))dt 

J  to  Jto  Jr\ 

=  9{tt)  [  sip(s(t  -  t0))dt 

Jto 


rto+S 

+g(t0  -f  £)  J  sip(s(t  -  t0))dt 
g{to)  /  ip(z)dz  +  g(t0  +  6)  f  ip(z)dz 
g(to)'<P+  +g(to  +  S) J  if)(z)dz 


=  g(to)il>+  as  s  -+  oo, 

where  we  have  used  the  definitions  of  ip-  and  ip+  given  in  Eq.  (22).  Similarly, 

fl° 

/  g(t)sip(s(t  —  t0))dt  — »  g(t0  )ip-  as  a  oo. 

Jto  -s 

Together  Eqs.  (29)  and  (30)  show  that 

y/sl2  -*  g{to  )ip-  +  g{t+  )ip+  <  oo  as  3—i  oo. 

This,  in  turn,  implies  that 

I2  — ►  0  as  s  — >  00. 

For  if  I2  — ►  K  ^  0  as  s  — ♦  00,  then  y/sli  — *  00  as  s  — ♦  00,  which  contradicts  Eq.  (31). 


(29) 

(30) 

(31) 

(32) 


Equations  (26)  and  (32)  prove  property  (i),  and  Eqs.  (28)  and  (31)  prove  property 
(ii).  If  g(t)  is  continuous  at  to,  then  g(t^ )  =  g(to  ),  and  property  (iii)  follows  immediately. 
This  completes  the  proof. 


Properties  (ii)  and  (iii)  imply  that  at  t  =  to  the  wavelet  transform  falls  off  to  0 
less  rapidly  along  the  s  suds  if  g(t)  has  a  jump  discontinuity  at  t0.  In  other  words,  any 
discontinuous  jump  of  the  function  g(t)  at  t  =  to  implies  the  value  of  <pg(s,  t0)  is  significant 
for  large  values  of  a.  In  a  practical  sense,  this  means  that  the  wavelet  transform  can  be 
used  for  transient  detection,  since  any  abrupt  change  or  short  term  feature  of  the  function 
should  cause  <pa(s,to)  to  be  significant  for  large  values  of  s.  Also,  because  of  the  assumed 
regularity  of  the  wavelet,  we  know  ^(0)  =  0,  hence 

0  =  $(0)  =  /  ip(t)dt  =  ip-  +  ip+  =>  ip+  =  —ip-.  (33) 

J  —  00 
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Furthermore,  the  condition  ip  6  Ll{ R )  implies  if  {  and  V’--  are  finite  since 

/°° 

hKOI<*  =  llV’lli  <  °o-  (34) 

-  oo 

Finally,  almost  any  function  encountered  in  practice  modeling  a  physical  phenomenon  will 
be  a  well  behaved,  bounded,  piece- wise  continuous  function.  Thus,  requiring  g{t)  to  be  of 
bounded  variation  is  not  overly  restrictive. 

Theorem  3:  Suppose  g(t)  and  the  wavelet  ip(t)  are  members  of  L2{R),  then 

l<M5>r)!2  <  Ml  HV'lll  (35) 


Proof  of  Theorem  3:  Equation  (35)  follows  directly  from  the  applica*:on  of  the 
Schwartz  inequality  to  Eq.  (5): 


|V>s(s,r)r  = 


< 


I  Vs  f  g(t)if*(s(t  -  T))dt 

I  J -OO 

[/oc  [S  ~  T^dt 

[/.«H  [/^  W(z)\2dz 

mi  mi 


(36) 


4.  LINEAR  SYSTEMS  AND  THE  WAVELET  TRANSFORM 

In  this  section  we  state  and  prove  three  results  of  the  wavelet  transform  as  applied 
to  linear  system  theory. 

Theorem  Let  x(t)  be  the  deterministic  input  to  a  linear  system  whose  known 
impulse  response  is  g(t),  and  whose  output  is  y[t )  =  x(t)  *  g(t),  where  *  denotes  the 
convolution  operator.  Furthermore,  let  the  wavelet  transforms  of  x(t),  g{t),  and  y(t)  be 
given  by 

W  jx(t)j  =  4>x{s,t)  -  y/s  J  x(t)i>*(s(t  T))dt, 

w  {</(*)}  -  M3’7)  ^  y/s  J  -  r))dt, 

W  |y(t)|  =  Ms’t)  ~  \T9  J  3/(0V’*(5(<  r))dt.  (37) 

Then  the  wavelet  transform  of  the  output  y(t)  is  related  to  the  wavelet  transforms  of  x(t) 
and  g(t)  through  the  equation 

/oo  fOC 

x(z)(j)g(s,z  r  )dz  =  /  g(z)(f)x(s,z  r)dz.  (38) 

■  OO  J  -  OC 
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Proof  Theorem  4:  First,  it  is  well  known  from  linear  system  theory  that 


y(t)  -  f  x(z)g{t  -  z)dz  -  x(t)  *  g(t).  (39) 

J  -  OO 

Substituting  Eq.  (39)  into  the  definition  of  the  wavelet  transform  of  y(t),  and  invoking  a 
change  of  variable  yields 

w  {y(o}  =  Msit) 

/OO 

y(t)ip*(s(t  -  r))dt 

-OO 

—  y/s  f  f  x(z)g{t  --  z)ij)*(s(t  T))dtdz 
J -OO  J  —  OO 

—  f  x(z)\s/s  f  g(t  z)rp'(s(t  -  r))dt  dz 

=  /  x(z )  /  y(0V’*(5(^  -  (T  -  z))dt  dz 

J -OO  J  —  oo 

=  /  x(z)4>g(s,T  -  z)dz,  (40) 

which  proves  the  first  integral  in  Eq.  (38).  The  second  integral  in  Eq.  (38)  is  derived 
similarly,  but  starts  with  the  alternate  form  of  the  convolution  integral 


/oo 

g{z)x{t  -  z)dz. 

-oo 


This  completes  the  proof. 


It  is  well  known  from  the  theory  of  stochastic  processes  that  the  autocorrelation 
function  of  a  stationary  process  x(t )  is  related  to  its  power  spectral  density  Sxx(f)  through 
the  Fourier  transform  relation 

E{x(t)x*(t  -  z)}  -  Rxx(z)  --  r  Sxx(f)e*"t*df.  (42) 

J  —OO 

This  is  known  as  the  Weiner- Khinchine  Theorem  [6].  Furthermore,  if  x(t)  is  passed  through 
a  linear  time-invariant  filter  whose  impulse  response  is  g(t),  then  the  power  spectral  density 
of  the  output  y(t)  [6]  is  given  by: 


2 

Syy(f)  =  T  {y(0} 


The  following  theorem  gives  similar  relationships  for  the  wavelet  transform. 

Theorem  5:  Let  x(£)  be  a  wide-sense  stationary  stochastic  process  whose  autocorre¬ 
lation  function  is  defined  by 

R~(z)  =  E{x(t)x'(t  z)\.  (44) 
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Then  the  wavelet  transform  of  the  autocorrelation  function  with  respect  to  the  wavelet 
■0(t)  is  given  by 

W  |iZx*(-2)}  =  4>xx(s,t)  =  Vs  J_  Rxx{z)i>'{s(z  -  r))dz,  (45) 

and  the  Fourier  transform  of  (/)xx(s,t)  with  respect  to  r  is 

T  {<M»,r)}  =  §„(*,/)  =  Wis^n,  (46) 

where  ^(/)  is  the  Fourier  transform  of  the  wavelet  Tp(t),  and  Sxx(f )  is  the  power  spectral 
density  of  ®(t).  Furthermore,  let  x(t)  drive  a  linear  time-invariant  system  whose  impulse 
response  is  g(t),  and  whose  output  is  denoted  by  y(t).  Then  the  Fourier  transform  of  the 
wavelet  transform  of  the  output  autocorrelation  function  Ryy(z)  is  given  by 

*vy(s,f)  =  \G(f)\2'S>xx(sJ).  (47) 


Proof  of  Theorem  5  (first  method):  We  take  the  wavelet  transform  of  the  autocor¬ 
relation  function  directly,  and  use  the  Fourier  transform  property  (a  transformation  with 
respect  to  the  variable  t) 


where  s  >  0.  Thus, 

W  |fixx(2)|  =  y/s  J  Rxx{z)^{s{t  -  r))dz 


=  s-1'2 


=  s 


=  a 


=  s"1/2 


r  rMT  y'(-f/*y2*f{z-T)df 

J—oo  |y — oo 

1/2  f°°  *•(_//,)[  r  Rxx(z)e**f*dz 

1/2  r  **(-//s)5xx(-/)e-'2^d/ 

7—oo 

J  —  OO 


dz 


2*fTdf 


Therefore, 


This,  in  turn,  implies  that 


which  proves  Eq.  (46). 


/  Sxx\ 

s/s 


*-(<•/)  =  y^-s„U), 


(48) 


(49) 

(50) 

(51) 
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If  we  now  consider  the  case  of  x(t)  driving  a  linear  system  whose  impulse  response 
is  then  from  Eq.  (43)  the  power  spectral  density  of  the  output  y(t )  is 


Syy{f)  =  \G(f)\2Sxx(f). 


Thus,  from  Eq.  (51)  and  (52)  it  follows  that 


*»(*,/) 


**(//*) 

*•(//-) 


Syy(f) 

\G{f)\*Sxx{f) 


s-tf) 


=  |G(/)|2f^^ 

Vs 

=  \G(f)\>*xx(s,f), 

which  proves  Eq.  (47).  This  completes  the  proof. 


Proof  of  Theorem  5  (second  method):  This  approach  uses  the  definition  of  the  auto¬ 
correlation  as  an  expectation.  Substituting  Eq.  (44)  into  Eq.  (45),  and  using  the  change 
of  variable  w  =  t  —  z  yields 

W  {#**(*)}  =  V»  f  E{x(t)x*(t  —  z)}tf>*(s(z  —  r))dz 
=  E^x(t)y/s  J  x*(t  -  z)ip*(s(z  -  r))dz| 

=  E|x(t)|\/ay  x(w)i>(s((t  —  t)  —  w))dw  | 

=  E{z(t)h*(t  -  t )},  (54) 

where  we  interpret  h(t)  to  be  the  response  to  the  input  x(t)  of  a  linear  time-invariant 
system  whose  impulse  response  is  v/i^(st).  Thus,  the  last  line  in  Eq.  (54)  is  the  cross- 
correlation  function  Rxh{T)-  The  Fourier  transform  of  this  correlation  function  is  the 
cross-power  spectral  density  function,  and  it  is  well  known  that 

T  (  D  /  —  'll  C  I  t\  ^  (/ /S)  o  /  t\ 


T  {iU(r)}  =  Sxh(f)  =  -^p-Sxx(f) 


where  we  have  used  the  Fourier  transform  relation 


|v/3V’(5t)| 


»(//«) 


We  note,  however,  that  Eq.  (55)  is  also  equal  to  the  Fourier  transform  with  respect  to  r 
of  the  wavelet  transform  4>xx{s,r ).  Thus,  we  have  proved  Eq.  (46).  The  proof  of  Eq.  (47) 
is  identical  to  the  approach  used  in  the  first  method.  This  completes  the  proof. 
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Let  R(s,t)  be  a  density  function  that  describes  the  average  amount  of  spread  in 
range  delay  r  and  Doppler  variable  s  a  signal  undergoes  when  applied  to  a  dispersive 
channel.  We  call  R{s,t)  the  channel  scattering  function  [7,8].  Furthermore,  if  the  input 
signal  to  the  channel  is  x(t),  and  its  wavelet  transform  is 

W  {x(t)>  =  (fx(s,T)  =  \fs  I  x(t)iJ>*(s(t-T))dt,  (57) 

J  —  OO 

then  the  expected  square  magnitude  of  the  wavelet  transform  of  the  channel  output  y(t) 
is  given  by 

E{\cf>y(s,T)\2}  =  [  f  R(w,z)Ax^(s/w,w(t  -  z))dwdz,  (58) 

J z=  —  oo  J w=0 


where  Ax^,(s,r)  is  the  wideband  crossambiguity  function  defined  by 


AXii,(s,T)  =  s 


[  x(t)ij)*(s(t  —  r))dt 

J  —  OO 


(59) 


Derivation:  In  a  channel  that  spreads  a  signal  in  both  time  and  Doppler,  the  output 
is  composed  of  time  shifted  and  Doppler  shifted  replicas  of  the  input  signal.  For  a  channel 
composed  of  a  finite  number  of  scatterers,  the  output  signal  would  be  given  by 

N 

y{t)  =  X  x(*(*  -  r0)>  (60) 

»= i 

where  pi  is  the  random  complex  reflection  coefficient  of  the  ith  scatterer.  By  Theorem 
1,  we  know  that  the  wavelet  transform  of  an  input  signal  dilated  by  an  amount  st-  and 
delayed  an  amount  T{  is  given  by 

W  {x(3i(f  -  T;))}  =  <t>x(s/Si,3i(t  -  rj).  (61) 

Thus,  by  linearity,  the  wavelet  transform  of  output  signal  in  Eq.  (60)  is 

w  {y(0}  =  <f>y{3iT)  =  '52pi<l>*{*/*ii*i{T  -  r«))-  (62) 

i= 1 

If  we  assume  that  the  expressions  in  Eqs.  (60)  and  (62)  are  approximations  for  a  channel 
described  by  a  continuous  distribution  of  scatterers,  then  as  N  — »  oo,  we  have 

[OO  [OO 

<f>y(s,r)  —  /  /  S(z,w)<j>x(s/w,  w(t  —  z))dwdz,  (63) 

J z— — oo  Jw=0 

where  5(s,r)  is  called  the  spreading  function,  and  describes  the  distribution  of  the  reflec¬ 
tion  coefficient  in  range  delay  r  and  Doppler  variable  s  [7].  It  is  a  stochastic  function. 
Furthermore,  assume  that  the  channel  exhibits  uncorrelated  spreading,  then 

E{5(s,r)5*(J,r)}  =  R(s,t)S(s  -  s)6(t  -  f),  (64) 
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where  £(•)  denotes  the  Dirac  delta  function.  We  apply  Eq.  (64)  to  the  correlation  function 
of  the  wavelet  transform  of  the  channel  output  as  follows: 


poo  p  oo  poo  poo 

E{<j>y(s,T)(f>l(s,T)}  =  I  J  j  I 

J Z—  —  oo  J w~0  *  z~  —  oo  Jw=0 

E{S(w,z)S*(w,z)} 

(f>x(s/w,w(T  -  z))<f>l(s/w,w(r  -  z )) 
duidwdzdz. 

poo  poo 
J  z—  —  oo  J w=0 

R(w,  z)(j>x(s/w,w(T  —  z))4>1(s/w,w(t  —  z))dwdz. 


(65) 


Setting  s  —  s  and  f  =  r  yields  Eq.  (58). 


5.  THE  WAVELET  EXPANSION  AND  DUABECHIES  WAVELETS 


The  transform  discussed  in  the  previous  sections  shows  that  for  a  suitable  choice  of 
V’(t),  any  function  in  L2(R )  can  be  expressed  as  an  integral  sum  of  time-dilated  and  time- 
delayed  wavelets.  Note  that  this  is  not  an  orthogonal  representation,  since,  in  general,  we 
do  not  have 

/<*> 

ip(ai(t  —  r1))-0*(s2(t  —  T2))dt  =  0,  for  all  ^  s2  and  rj  ^  r2.  (66) 

-OO 

Thus,  the  transformation  is  in  some  sense  redundant.  However,  an  orthogonal  represen¬ 
tation  can  be  derived,  provided  the  wavelet  ip(t)  obeys  some  additional  restrictions.  The 
result  is  an  expansion  composed  of  a  countable  sum  of  wavelets. 

The  description  of  this  expansion  was  first  given  by  Mallat  [5].  In  his  work,  he 
presented  the  concept  of  a  multiresolution  approximation  of  any  function  g(t)  £  L2(R). 
This  is  a  nested  sequence  of  closed  subspaces  in  L2(R)  denoted  as  { Vj}jez ,  where  Z  is  the 
set  of  all  integers  for  which  the  following  are  true: 


(a)  V,  C  Vj+U  for  all  j  £  Z, 

(b)  Ujez  vi  is  dense  in  L2{R)  and  f]jeZ  Vj  =  0, 

(c)  g(t )  £  Vj  g(2t)  £  Vj+1  for  all  j  £  Z 

(d)  g{t)  £  Vj  =>  g{t  —  2~*k)  £  Vj  for  all  j  £  Z 

(e)  Let  l2(Z)  denote  the  space  of  square  summable  infinite  sequences,  then  there  exists 

an  isomorphism  T  :  V0  — >  l2(Z)  which  commutes  with  the  action  of  Z. 
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Properties  (a)  and  (b)  merely  state  that  the  nested  sequence  of  subspaces  must  span 
the  space  L2(R).  Properties  (c)  and  (d)  describe  the  effects  of  time  dilation  and  time 
delay.  Collectively,  however,  all  the  properties  can  be  used  to  prove  that  there  exists  a 
function  0(£)  €  L2(R)  such  that  for  any  j  £  Z  the  subspace  Vj  is  spanned  by  the  set 
{s/2iO{V{t  -  k))}keZ  ■  0(t)  is  called  the  scaling  function ,  and  the  properties  stated  above 
imply  that  if  we  wish  to  approximate  (in  a  minimum  integral  square  error  sense)  the 
function  g(t )  by  a  function  gj(t)  in  the  subspace  V),  then 

&■(*)=  E  9i,kV226{V{t  -  r2k)),  (67) 

k—  —  oo 

where  the  coefficients  of  the  expansion  are  given  by 

gj<k  =  V2i  [°°  g(t)8\2j(t  -  T-2k))dt.  (68) 

This,  of  course,  is  reminiscent  of  a  generalized  Fourier  expansion;  however,  the  difference 
here  is  that  the  expansion  is  only  over  the  subspace  Vj  not  the  entire  space  L2(R).  Equation 
(67)  explains  property  (e)  in  that  the  sequence  {gj,k}ktz  is  an  alternate  representation  of 
g(t).  Because  the  set  {\/2i0(2 J(t  —  k))}kez  is  a  basis,  the  elements  are  orthogonal.  Thus, 
by  calculating  HgH’  one  finds  that  it  is  equal  to  the  sum  of  the  squares  of  the  sequence 
{ 9j,k}kez ■  Since  g(t)  e  L2(R),  ||g||£  is  finite  implying  {gj,k}kez  £  ^2(^)- 

For  practical  applications,  we  need  to  know  more  than  just  the  existence  of  the 
scaling  function;  we  need  to  know  how  it  is  parameterized,  and  how  to  compute  it.  A  step 
in  this  direction  is  the  following  theorem  that  is  proved  in  Refs.  2  and  5. 

Theorem  6:  Let  0(t)  be  a  scaling  function,  and  H(f)  be  the  Fourier  series  defined  by 

«(/)=  E  (69) 

k=—oo 

where  {/i(A:)}*ez  is  the  sequence  defined  by 

h(k)  =  -  r  8(2~1t)9*(t  +  k)dt.  (70) 

Z  J  —  OO 

Then  H(f)  satisfies  the  following  properties: 


(i)  |£T(0)|  =  1, 

(ii)  h(k )  ~  0(k2)  as  k  —*  oo, 
(m)\H(f)\2  f  |ff(/+ 1/2)|2  =  1. 


Furthermore,  let 


\H(f)\^0  for  /  e  [0,1/2), 


(71) 
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then  the  Fourier  transform  of  the  scaling  function  is  given  by 

€>(/)  =  n  H( 2  ”/)•  (72) 

P-1 

The  proof  exploits  properties  (a)  through  (e),  which  define  the  multiresolution  approxi¬ 
mation.  In  particular,  it  uses  properties  (c)  to  say  that  the  function  8(t/2)/2  is  equal  to  a 
weighted  sum  of  the  functions  9(t  +  A:),  and  the  weights  are  the  elements  of  the  sequence 
{/i(/c)}fcez.  Furthermore,  any  sequence  used  to  define  H(f)  in  Eq.  (69)  so  that  H(f)  obeys 
properties  (i)  through  (iii)  of  the  theorem  can  be  used  to  find  0(t)  by  defining  its  Fourier 
transform  through  Eq.  (72).  Thus,  our  choice  of  the  h(k)’s  is  somewhat  arbitrary. 

The  expansion  given  by  Eq.  (67)  is  simple,  but  not  always  convenient  for  practical 
applications.  For,  if  we  wanted  the  approximation  of  g(t)  in  the  subspace  V)  and  had  the 
expansion  for  the  approximation  in  V).  i,  we  would  have  to  recompute  all  expansion  coef¬ 
ficients.  Furthermore,  this  expansion  does  not  lend  itself  to  defining  filtering  operations. 
An  alternate  expansion  can  be  derived  by  noting  (through  the  projection  theorem)  that 
there  exists  a  subspace  Oj  composed  of  functions  that  are  orthogonal  to  those  composing 
Vj  such  that 

Oj®Vs  =  Vltl,  (73) 

where  ®  denotes  the  Cartesian  product.  Thus,  from  property  (b)  of  the  multiresolution 
expansion  one  can  show  that 

u°i  =  (74) 

)62 

In  light  of  this  new  definition  we  have  the  following  theorem  whose  proof  can  be  found  in 
the  Refs.  2  and  5. 

Theorem  7:  Let  {Vj}jeZ  define  the  multiresolution  approximation  of  the  space  L2(R), 
6(t)  be  the  scaling  function  whose  Fourier  transform  is  ©(/),  and  H(f)  be  the  Fourier 
series  describing  the  Fourier  transform  of  9(t)  as  in  Theorem  6.  Then  there  exists  a 
function  i>(t)  such  that  {'/2Jtp(2 i(t  -  2 ~2k))}j<k^ZyZ  is  a  basis  for  L2(R),  and  the  Fourier 
transform  of  i^(t)  is  given  by 


9{f)  -  k(£)q( 

\  2  /  V  2  / 

(75) 

where 

K(f)  =  e  ]2lrf H*(f  f  1/2). 

(76) 

From  Eqs.  (75)  and  (76)  it  is  possible  to  show  that  6(t)  is  equal  to  a  linear  sum  of  time 
delayed  scaling  functions.  The  function  is  called  an  orthogonal  wavelet,  and  the 
theorem  given  above  says  that  any  function  g(t)  in  L2(R)  can  be  written  as 


<7(0  =■  £  <7j,fc  v^(2j(t  -  2  -*'*)),  (77) 

j,kf  7j  x  Z 

where 

9 fX  g{tW(2’{t  2  >k))dt  =  <f>g(22 ,2~}k).  (78) 
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Thus,  we  see  that  the  wavelet  expansion  coefficients  are  equal  to  the  value  of  the  wavelet 
transform  at  the  points  (s,r)  =  (2J,  2_J/c).  Figure  1  shows  this.  Moreover,  we  see  that 
6[t)  plays  a  role  in  the  definition  of  "0(0-  From  a  computational  standpoint,  Theorem  7 
says  that  we  must  find  6{t)  (or  its  Fourier  transform)  first,  and  then  compute 


0  12  3  4 


X 

Fig.  1  —  Points  in  the  (.s.  r)  plane  where  the  wavelet  expansion 
coefficients  are  equal  to  the  wavelet  transform.  This  figure  shows 
all  points  lying  within  and  on  the  boundary  of  the  region 
s  €  [2~5.  2:|  and  r  =  10.  4], 


Theorems  6  and  7  not  only  suggest  how  the  scaling  function  and  orthogonal  wavelets 
can  be  computed,  but  also  suggest  that  we  can,  to  some  degree,  control  the  shape  of  the 
wavelet  in  the  time  domain  according  to  how  we  choose  the  sequence  {h(fc)}*e£.  One 
desirable  property  is  to  have  a  wavelet  with  compact  support  in  the  time  domain,  i.e., 
it  is  time  limited  in  that  it  is  nonzero  only  over  a  given  interval.  Such  a  wavelet  gives 
a  true  sense  of  time  locality.  A  set  of  orthogonal  wavelets  with  compact  support  was 
discovered  by  Daubechies  [2].  They  are  parameterized  by  an  integer  n,  are  real  valued, 
and  are  denoted  as  Vv.(0  for  n  >  2.  In  fact, 


aupp  ifin  C  1(1  -»),»], 


f  y>„(2’(t  —  2  */)-0^(2,7(t  —  2  *k))dt  =  0  for  all  /  ^  k  or  i  /  j. 

J -oo 
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The  associated  scaling  function  is  denoted  by  6n(t).  These  wavelets  are  derived  by  choosing 
the  sequence  {h(k)}keZ  so  that  is  of  finite  length  (a  FIR  filter  in  signal  processing  parlance). 
The  result  is  the  set  of  sequences  {hn(k)}k£Z  that  are  nonzero  for  fc  =  0, 1, . . . ,  (2n  —  1); 
thus,  they  are  of  length  2n.  Details  of  the  procedure  for  finding  these  sequences  can 
be  found  in  Daubechies’  original  paper;  however,  we  have  tabulated  the  sequences  for 
n  =  2, 3, . . . ,  14  in  the  appendix. 


Figure  2  shows  the  Daubechies  orthogonal  wavelets  for  n  =  2, 3,...,  13  They  were 
generated  by  first  calculating  (approximating)  the  Fourier  transform  of  the  associated 
scaling  function  via  the  equation 


where 


=  II  ffn(2 

p=0 

(81) 

2n~l 

£  hn{k)e-*'kl. 

(82) 

k=o 


Once  ©n(/)  is  found,  we  calculate  the  Fourier  transform  of  the  orthogonal  wavelet  as 


*n(/)  =  ^({)©„(0,  (83) 

where 

*»(/)  =  +  1/2).  (84) 

Equations  (81-84)  follow  directly  from  Eq.  (69),  (72),  (75),  and  (76).  In  particular,  the 
truncated  product  in  Eq.  (81)  gives  good  results  for  P  —  20  for  low  values  of  n  (n  =  3), 
to  P  =  25  for  high  values  of  n  (n  =  13).  This  was  checked  by  calculating  the  normalized 
cross  correlation  between  two  Daubechies  wavelets  of  order  n,  where  one  was  derived  by 
using  P  =  N,  and  the  other  with  P  =  N  +  1.  For  P  =  25  (or  P  =  20  for  low  values  of  n) 
the  correlation  was  negligibly  different  from  1. 


6.  PROPERTIES  OF  DAUBECHIES  WAVELETS 


We  now  state  and  prove  five  theorems  about  the  Daubechies  wavelets  introduced 
in  the  previous  section.  The  first  three  theorems  state  that  these  wavelets  are  bounded, 
continuous,  and  in  most  cases  differentiable. 

Theorem  8:  ^>„(£)  is  bounded  for  all  n. 
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8  -6  4  2  0  2  4  6  8  -8  6  -4  -2  0  2 

TIME  I  TIME  I 


Fig.  2b  —  Daubechies  orthogonal  wavelets  for  n  =  8,  . . .  ,13.  These  wavelets  are 
nonzero  only  in  the  interval  [(I  -  n),  n]. 
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Proof  of  Theorem  8:  Daubechies  showed  that  the  Fourier  transform  of  the  wavelet 
ipn(t)  has  the  property 


|©„(/)i  <  c(  1  +  |/|)-"+ios£/i°s2  <  c  <  oo, 

for  some  constant  C,  and  where  B  <  2n_1.  It  follows  that 

l«»(OI  =  I  r  Qn(f)e’2*ftdf 

\J  —  OO 


(85) 


-oo 

r  oo 


<  c 


r  \On(f)\df 

J  —  oo 

f  df 


(1  _|_  \f^n-logB/log2  ’ 


and 


log  B 


<  n  —  1  =>  r*  - 


log  B 


>  1. 


log  2  log  2 

Therefore,  there  exists  some  e  >  0  such  that  Eq.  (86)  can  be  rewritten  to  yield 

df 


(86) 


(87) 


l«»MI  <  C  f 


=  2C 


=  2  C 


OO 

OO 


(1  +  l/l)1+e 

Jf _ 

o  (1  +  /)1+£ 

dz 


rl+c 


=  -2  C- 


2  C 


<  oo. 


(88) 


Thus,  the  scaling  function  is  bounded  for  all  t.  Since  {hn(k)}kez  is  a  finite  sequence, 
it  is  possible  to  show  through  Eqs.  (81)  and  (82)  that  ^„(t)  is  a  finite  sum  of  time  delayed 
scaling  functions.  Therefore,  the  wavelet  0n{t)  is  bounded  for  all  t. 


Theorem  9:  ^„(t)  is  continuous  for  all  n. 


Proof  of  Theorem  9:  Let  Ca(R )  be  the  set  of  all  functions  such  that 

g(t)  £  Ca(R)  «•  /~  |G(/)|(1  +  |/|)*+“  <  oo.  (89) 

If  a  =  k  for  k  =  0, 1,2, . . .,  then  Cfc(i?)  is  the  space  of  k-times  continuously  differentiable 
functions,  where,  in  particular,  C°(R )  is  the  space  of  continuous  functions.  Daubechies 
has  already  shown  that  ipn{l)  £  Ca(R)  for  some  a  >  0.  Thus,  we  only  need  to  establish 
that  this  implies  ip(t)  (E  C°(R).  To  begin,  note  that  for  any  e  £  [0,  a] , 


20 


NRL  REPORT  9316 

(1  +  l/l)1+<  <  (1  +  l/l)1+“,  for  all  /.  (90) 

This  implies  that 

IG(/)|(1  +  l/l)1+‘  <  |G(/)|(1  +  |/|)1+“,  for  all  /,  (91) 

which,  in  turn,  implies 

r  |G(/)|(1  +  |/|)1+*<i/  <  r  |G(/)|(1  +  \f\y+°df  <  00.  (92) 

J~oo  J— oo 

Thus,  by  the  definition  of  Ca(R)  as  it  follows  from  Eq.  (89),  we  see  that  <£„(£)  £  Ce(R) 
for  all  e  €  [0,aj.  Therefore,  V’n(t)  £  C'°(iE),  and  so  ipn(t)  is  continuous. 

Theorem,  10:  ipn(t)  is  continuously  differentiable  for  n  >  4. 

Proof  of  Theorem  10:  Daubechies  proved  that  a  >  1  for  n  >  4.  Furthermore,  from 
the  proof  of  the  previous  theorem,  we  know  that  £  Ce  for  all  positive  e  less  than  a. 

Thus,  V’n(t)  £  C1(R)  for  all  ra  >  4.  In  other  words,  ifn(t)  is  continuously  differentiable  for 
n  >  4. 

Theorem  11:  For  n  >  4,  if>n{t)  has  a  finite  spectral  variance,  i.e., 

/OO 

/2|*n(/)|d/<oo.  (93) 

-OO 


Proof  of  Theorem  11:  For  n  ^  4,  rfn(t)  is  continuously  differentiable,  therefore, 
is  continuous.  Also,  since  has  compact  support  (an  interval),  so  does  Tp'n(t). 

Furthermore,  since  a  continuous  function  over  a  compact  interval  is  bounded,  we  see  that 
ip'n(t)  is  bounded  over  supp  ipn,  thus 


r  m^dt  =  [  wn(t)\2dt 

J  —  oo  J  »utm  t br%. 


,upp  \/>n 

<  meas{supp  ifn)  •  max  \ip'n(t)\2 

lupp 

<  meos([(l  —  7i),7i])  ■  max  \ij>'n{i)\2 

te[(l-n),n] 

=  (2n  —  1)  •  mai  |0)iJ  <  oo, 

te[(l-n),n] 


(94) 


where  meas(-)  is  the  Lebesgue  measure,  and  have  used  Eq.  (79).  Since  the  Fourier  trans¬ 
form  of  Tp'n(t)  is  j2t /$„(/),  we  have  by  Parseval’s  Theorem, 


4tt2  rf2\nf)\2df=  r  mt)?dt 

J  —  oo  J -oo 


<  oo. 


(95) 
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Theorem  12:  For  n  >  4,  or  for  any  wavelet  with  a  finite  spectral  variance,  if  g(t)  £ 
L2(R),  then 

£M£iI)  '  _  0,  as  .  0.  (96) 


Proof  of  Theorem  12:  By  the  definition  of  the  waveL  transform  we  haw 

^  -  r))it 

=  - s W  P  g(t)rMt  -  r))dt.  (97) 

J  ~oo 

Therefore,  taking  the  square  magnitude  of  Eq.  (97),  applying  the  Schwartz  inequality,  and 
using  Eq.  (95)  from  Theorem  11,  we  have 

<  *  r  \g{t)\2itP  -T))fit 

J  —  oo  J  —  oo 

=  ,2  r  \g{t)\2dt  r  \rM\2dz 

J  —  OO  j  —  OO 

/OO 

/2|*n(/)|2d/<oo.  (98) 

-oo 

Thus,  by  applying  the  squeeze  theorem  for  limits,  Eq.  (96)  follows  from  Eq.  (98).  This 
proves  the  theorem. 


Basically,  this  theorem  says  that  as  s  approaches  zero,  the  derivative  of  the  wavelet 
transform  along  the  r  axis  approaches  zero.  In  other  words,  <P(s,t)  becomes  ‘smoother’ 
along  r  as  s  — >  0. 

Theorem  13:  ^„(t)  £  i1(/Z)  fj  L2(R)  for  all  n. 


Proof  of  Theorem  13:  By  Theorem  8  we  know  that  ifn{t)  is  bounded,  i.e.,  |^>n(i)|  <  K 
for  some  K  <  oo.  Therefore,  using  Eq.  (79)  we  have 

f  \r[>n(t)\dt  <  meas  ( supp  xfn)  •  K  <  (2 n  —  1)  •  K  <  oo,  (99) 

J  —  OO 


which  implies  that  ipn(t)  £  T1(fi).  Also, 


/  \i>Ti{t)\2 dt  <  meas  ( supp  tfn)  •  K 2  <  (2n  —  1)  ■  K2  <  oo, 

J  —  OO 


(100) 


which  implies  that  i>n(t)  £  L2{R).  This  proves  the  theorem. 
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This  last  theorem  shows  that  the  Daubechies  orthogonal  wavelets  obey  the  regularity 
property  as  defined  in  Eq.  (9).  This  says  that  the  wavelet  transform  that  uses  a  Daubechies 
orthogonal  wavelet  as  the  transform  kernel  is  meaningful  in  the  sense  that  its  inverse  exists. 


Proof  of  Theorem  1 4:  Daubechies  and  Mallat  have  shown  that 


OO 


*»(/)  =  +  1/2)  n  H.UI 2‘). 

fc=l 

(101) 

where 

»„(/)=  + 

(102) 

and 

=  E  (n  ’  ‘ +  *)  »»“(*/)• 

(103) 

From  Eqs.  (102)  and  (103)  one  finds 

/im|tfn(/)|2  =  1, 

(104) 

which  implies 

OO  2 

lim  n  =1. 

'  fc=l 

(105) 

Now  consider  the  limit 


/— *o  I/I  /-o 


(1  -  e~^f) 
2 


IQn(/+l/2)|2 


(106) 


Using  Eq.  (103),  and  L’Hopital’s  rule,  we  see  that 

Um  '9-(/+ 1/2)!!  .  lim 

/-o-  |/|  /->0+  Djf 

=  2 

h  v  * 

=  0. 


sin2k  1('n'f  +  7t/2)  cos(irf  -f  7t/2) 


(107) 
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Similarly, 


Um  l QM+llDL  _  o. 
/-«-  i/l 


(108) 


Combining  Eqs.  (101-108)  yields 


I^(/)|2  _  |^n(/+l/2)|2  .ft TirflokS 

lji™  i/i  ~  lf  ™  i/i  '^nw/2) 

J(l-e-^/)|2n  |Q„(/+  1/2)|2 


=  lim  — - - - L  •  lim  ^  w  - 

/-°  2  /-o  |/| 

=  001  =  0, 


■  n  j«(//2*) 
~ 1 _ 1 


(109) 


which  proves  (i).  (ii)  follows  from  (i)  for,  if  |^„(/)|2  — »  if  ^  0  as  /  — ►  0,  then  the  limit  in 
Eq.  (109)  would  be  infinite. 


(iii)  is  proved  in  two  parts,  since  we  can  write  the  integral  in  (iii)  as 

r\MniiJ  =  f  +  r  M,/  =  a  +  /2. 

J-oo  |/|  J- 1  J\f\>l  |/| 

Consider  I\.  We  first  note  that  because  |.ffn(/)|2  <  1,  and  using  Eq.  (101),  that 


(110) 


l*n(/)|J  =  |e-'!’/|2|H„(/+l/2)|Jn^(//2‘)  <|ff„(/+  1/2)|!.  (Ill) 


Also,  from  Eqs.  (102)  and  (103)  we  have 


1  4-  e~j2irf  2n 

\Hn(f)\ 2  =  -^1 -  |C?n(/)|2 

1  _L  eJ2*7  2" 

5  -V- 

=  «*’>/)  ■  |<?„(/)|J. 


(112) 


Therefore,  since  sin(x)  <  |x|,  from  Eq.  (112)  it  follows  that 


l^n(/  +  l/2)|2 

I/I 


<  max  \Qn(f)\2 
/€[-  1/2.1/2J 


max  |C?n(/)|: 
/e[— i/a,i/2] 


<  max  |Q„(/)|‘ 

/e[- i/2,i/2]  ^ 


21  cos2n(7r/  +  tt/2) 

I/I 

2  sin2n(irf ) 

J  I/I 
,1  (^/)2W 
.  I/I 


(113) 
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Integrating  both  sides  of  Eq.  (113)  yields 


i  | Hn(f  +  1/2) j 


,ri";<00'  (114) 


which  proves  that  I\  <  oo.  Note  that  for  |/|  >  1, 


!^n(/)|2 


<  l*n(/)l2 


(115) 


Since  V’n(i)  £  L2(I2)i  by  Parseval’s  Theorem  we  know  that  $„(/)  €  L2(R),  so  it  follows 
that 

h = /  <  r  i»»(/)ia4f  <  oo.  (He) 

/m>i  /  J-  oo 


/|/l>i  I/I 

Combining  Eqs.  (114)  and  (116)  yields 

r  !*«(/)!’ 


d/  =  Ix  +  I2  <  OO, 


(117) 


which  proves  (iii). 


7.  APPLICATION:  WAVELET  PASSBAND  FILTERING 


Frequently  one  must  filter  a  passband  signal  to  reduce  transients  or  reject  out  of 
band  signals.  Generally,  this  is  accomplished  by  passing  the  signal  through  an  analog 
filter  designed  to  pass  only  those  frequency  (Fourier)  components  that  occupy  the  signal 
passband.  This  is  equivalent  to  convolving  the  input  signal  with  the  impulse  response  of 
the  filter.  In  this  section,  we  show  how  filtering  can  also  be  accomplished  by  using  the 
wavelet  expansion. 

It  was  shown  in  Section  5  that  if  g(t)  £  L2(R )  has  the  wavelet  expansion 

g(t)a=  £  gjtk  V2^{V{t  -  2-^fc)),  (118) 

jtk£Z  x  Z 

where 

gjtk  =  y/&  g(t)r{V{t  ~  2 ~3k))dt.  (119) 

We  realize,  however,  that  the  components  of  the  expansion  associated  with  large  values 
of  2-*  roughly  correspond  to  short  term  features  (transients)  of  the  signal  g(t),  and  the 
components  associated  with  small  values  of  2^  contribute  to  the  long  term  (average  or  dc) 
components  of  g(t).  Thus,  to  do  passband  filtering  in  the  wavelet  domain,  we  can  construct 
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a  signal  with  a  wavelet  expansion  that  uses  only  those  coefficients  that  correspond  to  a 
narrow  range  of  dilations.  In  mathematical  terms,  we  construct  the  signal 

MO  -EE  9:,k  VHil>{2  j(t  -  2 -ik)).  (120) 

j~jl  k  £Z 

Thus,  we  are  constructing  a  signal  with  the  coefficients  associated  with  the  wavelets  whose 
dilations  are  2J| , . . . ,  2>K . 


We  demonstrate  this  filtering  method  by  the  example  that  follows.  Figure  3  shows 
a  real  valued  signal  composed  of  two  windowed  sine  waves  and  a  spike  (transient).  The 
uncontaminated  signal  is  mathematically  expressed  as 

g(t)  =  i  exp((t  -  t)2/(rlj  [ sin(2irfit )  +  sin(27r/2t)],  (121) 


where  t  =  8.533,  cr i  =  3.413,  fi  =  1.5  Hz,  and  /2  =  2 f\  —  3.0  Hz.  Whereas  the  spike  is 
given  by 

n(t)  =  3exp(^(t  -  i)/<r£j,  (122) 

where  cr2  =  0.071.  Figure  4  shows  the  Fourier  spectrum  of  g(t )  +  n(t),  and  Fig.  5  shows 
its  wavelet  transform.  The  vertical  ridge  in  Fig.  5  is  due  to  the  transient.  Clearly,  the 
wavelet  transform  show  the  locality  of  the  transient  in  time,  which  the  Fourier  spectrum 
does  not.  We  first  analog  filter  the  signal  in  Fig.  3  by  passing  it  through  a  Chebyshev 
passband  filter  whose  Fourier  spectrum  is  given  by 


C(f)  = 


1 

S4  +  1.8037753  +  2.62680S2  +  2.025505  +  0.82851  ’ 


(123) 


where 

p  -  P 

S  =  ]2lZ  Bf'  (124) 

and  f0  =  \/ fhfi  is  the  geometric  mean  of  the  filter  passband  whose  width  is 

B  =  fh  —  fi  where  fi  <  fh-  In  this  case,  /;  =  1.15  Hz,  and  //,  =  3.35  Hz.  The  Fourier 
spectrum  of  the  filter,  shown  in  Fig.  6,  displays  some  ripple  in  the  passband.  This  is  an 
inherent  characteristic  of  the  Chebyshev  filter  type  [9],  For  this  particular  filter,  the  ripple 
width  (peak-to-peak  difference)  is  0.1  dB.  Figure  7  shows  the  Chebyshev  filter  output  in 
time  in  response  to  the  signal  in  Fig.  3,  and  Fig.  8  shows  its  Fourier  spectrum.  In  Fig.  7, 
we  note  two  features.  First,  we  see  that  as  compared  to  the  original  function,  the  two 
sine  waves  are  displaced  in  time.  This  is  due  to  the  phase  characteristic  of  the  filter: 
each  sinewave  experiences  a  different  phase  shift.  Second,  we  see  that  the  entire  signal  is 
displaced  in  time  slightly,  thus  accounting  for  a  group  delay  that  is  also  due  to  the  phase 
characteristic  of  the  filter. 


TIME  I 


FREQUENCY  I 


3  —  Real  valued  input  signal  composed  of 
two  windowed  sine  waves  and  a  spike 


Fig.  4  —  Square  magnitude  of  the  input 
signal’s  Fourier  spectrum 


Fig  5  --  Square  magnitude  of  the  input  signal's  wavelet  transform. 
The  transform  kernel  is  a  Daubechies  orthogonal  wavelet  of  order 
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FREQUENCY  * 

Fig.  6  —  Square  magnitude  of  the  Chebyshev 
passband  Filter's  Fourier  spectrum 
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TIME  I 

Fig.  7  —  Output  of  the  Chebyshev  passband  filter  in 
time  in  response  to  the  signal  in  Fig.  3 


Fig.  8  —  Square  magnitude  of  the  Fourier  spectrum  of 
the  input  signal  after  analog  filtering  with  a  Chebyshev 
passband  filter 


Figures  9  and  10  show  the  result  of  applying  the  wavelet  filtering  method  to  the 
signal.  Figure  9  shows  the  result  for  a  reconstruction  using  dilations  21,  22,  and  23.  Figure 
10  shows  a  reconstruction  using  slightly  less  wavelet  bandwidth.  In  this  case,  we  have 
used  dilations  21  and  22.  As  compared  to  passing  the  signal  through  a  Chebyshev  filter, 
we  see  qualitatively  that  wavelet  filter  produced  more  amplitude  distortion  in  the  time 
domain  but  does  not  show  the  same  effect  of  phase  distortion  or  group  delay.  Both  analog 
filtering  (using  the  Chebyshev  filter)  and  wavelet  filtering  de-emphasize  the  spike  but  do 
not  completely  remove  it. 

Figures  11  and  12  show  the  Fourier  spectrum  of  the  signals  derived  through  wavelet 
filtering.  They  show  that  the  low  frequency  components  of  the  spectrum  have  been  re¬ 
moved,  and  this  is  reasonable  since  a  large  part  of  the  spectral  energy  of  the  spike  is 
located  there.  Furthermore,  we  see  that  in  Fig.  12  a  spurious  peak  occurs,  which  implies 
that  wavelet  filtering  causes  nonlinear  distortion  in  the  Fourier  frequency  domain.  Figures 
13  to  15  show  the  wavelet  transforms  of  the  signals  resulting  from  wavelet  filtering,  and 
by  analog  filtering.  (Remember  that  the  wavelet  expansion  coefficients  of  these  functions 
are  equal  to  their  wavelet  transforms  at  the  points  (s,r)  =  (2%  2 ~lk)  for  t,  k  £  Z  x  Z .) 
As  compared  to  the  wavelet  transform  of  the  input  signals  shown  in  Fig.  5,  these  figures 
show  a  reduction  of  the  ridge  in  the  (s,r)  plane  associated  with  the  spike. 
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Fig.  9  —  Wavelet  reconstruction  of  the  input  signal  using 
dilations  2'  2:.  and  2\  The  wavelet  expansion  used  the 
Daubechies  orthogonal  wavelet  of  order  n  =  8. 
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Fig.  10  —  Wavelet  reconstruction  of  the  input  signal  by 
using  dilations  21  and  22 .  The  wavelet  expansion  used 
the  Daubechies  orthogonal  wavelet  of  order  n  —  8. 
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Fig.  1 1  —  Fourier  spectrum  (square  magnitude)  of  the 
wavelet  reconstruction  of  the  input  signal  using  dilations 
21,  2:,  and  2'.  The  wavelet  expansion  used  the 
Daubechies  orthogonal  wavelet  of  order  n  =  8. 


Fig.  12  —  Fourier  spectrum  (square  magnitude)  of  the 
wavelet  reconstruction  of  the  input  signal  using  dilations 
21  and  2:.  The  wavelet  expansion  used  the  Daubechies 
orthogonal  wavelet  of  order  n  =  8. 
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Fig.  13  —  Wavelet  transform  (square  magnitude)  of  the  wavelet 
reconstruction  of  the  input  signal  using  dilations  21,  22,  and  23.  The 
transform  kernel  is  a  Daubechies  orthogonal  wavelet  of  order  n  =  8. 


Fig  14  —  Wavelet  transform  (square  magnitude)  of  the  wavelet 
reconstruction  of  the  input  signal  using  dilations  2'  and  22.  The 
transform  kernel  is  a  Daubechies  orthogonal  wavelet  of  order  n  -  8 
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Fig  15  —  Square  magnitude  of  the  wavelet  transform  of  the  input 
signal  after  analog  filtering  with  a  Chebyshev  passband  filter.  The 
transform  kernel  is  a  Daubechies  orthogonal  wavelet  of  order  n  =  8. 


8.  APPLICATION:  A  DECONVOLUTION  ALGORITHM 


The  deconvolution  of  two  signals  has  important  applications  in  geophysics  and  com¬ 
munication  systems.  The  problem  can  be  simply  stated  as  follows:  given  a  known  input 
signal  x(t )  driving  a  linear  system  whose  unknown  impulse  response  is  #(£),  estimate  ^(t) 
given  the  known  (measured)  output  signal  y(t)  =  x(t)  *  g(t). 


The  easiest  way  to  deconvolve  is  to  simply  calculate  the  Fourier  transforms  X(f) 
and  Y{f),  and  find  their  quotient  G(/)  =  Y(f)fX{f).  This  is  reasonable  in  principle  but 
in  practice  can  produce  numerically  unstable  results,  since  one  may  have  to  divide  Y(  f)  by 
a  very  small  X(f)  if  our  discrete  numerical  approximation  of  X(f)  brings  us  close  to  one 
of  its  zeros.  This  problem  may  become  worse  in  the  presence  of  noise  in  the  measurement 
of  X{f).  The  following  outlines  a  deconvolution  procedure  in  the  context  of  the  wavelet 
expansion. 


by 


Both  the  impulse  response  g{t)  and  the  output  y(t)  have  wavelet  expansions  given 


g(t)  =  £  ***(*(*  -  »2 -■)), 

i,n£Zx  Z 

y(t)  =  £  m2  ■')), 

ji,m€  Zx  Z 


(125) 
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where  as  compared  to  the  definition  in  Eq.  (77)  we  have  dropped  the  factors  \/2*  and  \/2^ 
inside  the  summations  for  notational  convenience,  i.e.,  s/2x  — >  and 

yjtm  V#  — ♦  yj,m-  The  coefficients  gi>n  are  unknown;  they  are  the  ones  we  wish  to  es¬ 
timate.  On  the  other  hand,  the  coefficients  yJim  are  known  and  are  calculated  by  the 
direct  application  of  Eq.  (78)  on  the  known  (measured)  output  y(t ).  We  also  know  that 

y(t)  =  x(t)*g(t) 

=  x(t)*  52  yi,nV»(2‘(f  -  n2“‘)) 

t,7 i£Z  x  Z 

=  52  9i,nx(t)*ii>{2i(t-n2~i)).  (126) 

t,n6  Z  x  Z 

Each  term  inside  the  last  summation  in  Eq.  (126)  has  its  own  wavelet  expansion  given  by 

x(t)*^(r(t-n2^))  =  52  0(2^  -  m2~i)).  (127) 

JjTTl £.Z  X  Z 

By  substituting  Eq.  (127)  into  Eq.  (126),  and  with  rearrangement  of  the  summations  we 
find  that 

*(<)*»(<)=  E  E  4m  *.»  V>(2'(!  -  m2-')).  (128) 

j,m£ZxZli,n€ZxZ 

Comparing  Eq.  (128)  with  the  second  line  of  Eq.  (125)  shows  that 

Vj,m  —  52  9i,ni 

t,n£ZxZ 

9i,-l 
9i,  0 
9i,  1 

This  suggests  the  following  expression: 

Y  =V  G,  (130) 

where  V  is  a  tensor,  and  Y  and  G  are  matrices.  Assuming  a  suitably  defined  inverse  of  T> 
exists,  the  wavelet  coefficients  for  the  expansion  of  g(t)  can  be  found  from  the  expression 

G=V~lY.  (131) 

A  time  series  for  the  impulse  response  g(t)  can  be  found  immediately  by  using  coefficients 
g,in  in  Eq.  (77). 

At  this  point,  some  comments  are  in  order.  First,  the  matrices  and  tensors  in 
Eq.  (130)  contain  an  infinite  number  of  elements.  However,  in  practice  one  should  find 
that  elements  will  approach  zero  as  the  magnitude  of  their  indices  i,  j,  m  and  n  become 
large.  This  follows  from  Theorem  2,  which  says  the  wavelet  transform  decays  as  the 
magnitude  of  s  and  r  become  large.  Since  the  coefficients  of  the  wavelet  expansion  are 
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equal  to  the  wavelet  transform  at  each  point  where  (s,r)  =  (2ji,  2 ~*k),  we  know  that 
the  coefficients  will  decay  to  zero  as  the  magnitude  of  j  and  k  become  large.  Thus,  for 
practical  purposes,  we  need  only  use  matrices  and  tensors  with  finite  numbers  of  significant 
coefficients. 

Given  that  we  will  truncate  the  tensors  and  matrices  defined  above,  the  problem  can 
be  recast  as  a  standard  linear  programming  problem  using  matrices  and  vectors.  First, 
from  Theorem  2  we  know  that  the  coefficients  for  the  wavelet  expansion  of  the  input  x(t) 
will  only  be  significant  for  limited  ranges  of  i  and  n ,  i.e.,  i  =  i/, . . .  ,i/,,  and  it  —  n;, .  .  .  ,  n/,. 
From  Eq.  (127)  and  the  orthogonality  of  the  wavelets  in  the  expansion,  we  know  that  we 
need  only  consider  those  elements  of  V  associated  with  these  values  of  ?  and  n.  Similarly, 
the  coefficients  for  the  measured  output  will  only  be  significant  for  a  range  of  j  and  m, 
i.e.,  j  =  ji,...,jh ,  and  m  =  (More  generally,  (n/,n/,)  and  (m^m/,)  could 

depend  on  i  and  j  respectively.  This  was  done  in  the  numerical  example  that  follows.) 
From  Eq.  (125),  it  is  implied  that  we  need  only  consider  those  elements  of  V  that  are 
related  to  the  coefficients  for  the  output  associated  with  these  values  of  j  and  m.  With 
these  facts  in  mind  and  by  using  Eq.  (130),  this  implies  that  the  problem  of  finding  the 
wavelet  expansion  of  the  system  impulse  response  (finding  the  coefficients  p,,n)  can  be  well 
approximated  by  the  vector  equation 

y  =  Dg,  (132) 

where 

<T‘Z  ...  ...  C'Zh 

Jli™!  Ji,rni  Jl.mj 


/Ph'n'  rP’K’nh 

jl,mh  •  •  •  uji,mh  ’  ’  ’  Uji,mh  •  •  •  U},,mh 


J*t  »nl  A'l'Tlh  rPh,Til  jih>nh 

•  •  •  ajh<mi  ’  •  ■  ujh,mi  •  •  •  ujh,m, 


fPl,nh  J ‘h,n‘  if’1 l'nh 

Jh  ••• 

V  =  (  Vjtxmi  •  ■  ■  Vii^rn^  ■  ■  •  Jbh.mi  •  •  •  l/j^TOn  )  i 

9  ~  (  9ii,ni  •••  9ik,m  9*h,nh)-  (13«1) 

We  see  that  V  is  a  matrix  of  size  (jh  -  jt  +  l)(mh  -  m/  +  1)  x  (tj,  —  it  +  l)(n/,  -  nt  +  1),  and 

the  vectors  y  and  g  are  of  length  ( jh  ~  ji  +  l)(^n6  —  m;  +  1)  and  ( ih  -  it  +  1  )(n^  —  n/  +  1) 

respectively. 

If  N  =  (jh  —  ji  +  l)(mh  -  mi  +  1)  =  (ih  -  it  +  l)(nh  —  nj  -f  1),  then  D  is  a  square 
matrix,  and  the  vectors  are  now  elements  of  the  Euclidean  space  RN .  We  also  see  that 
the  solution  g  does  not  change  if  we  multiply  both  sides  by  the  Hermitian  of  D,  since 

6  =  DHy  =  DHDg  =  Qg.  (134) 
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We  now  seek  the  solution  to  the  equation 


6  =  Qg, 


(135) 


where  Q  =  DH D  is  now  a  positive  definite  matrix.  It  is  now  possible  to  show  that  finding 
the  solution  to  Eq.  (135)  is  equivalent  to  finding  the  solution  of  the  quadratic  problem 


min 

geR" 


\aHQa  -  f>"s 


(136) 


Finding  g  is  now  recast  as  a  minimization  problem,  thus  allowing  the  use  of  any  method 
designed  to  solve  this  class  of  problem. 


We  could  find  g  by  multiplying  b  by  the  inverse  of  Q,  providing  Q  is  nonsingular.  If 
it  is,  then  Q_1  is  unique,  and  we  have  a  unique  solution  for  g.  However,  rather  than  finding 
Q~l  directly,  we  really  need  only  to  find  g  directly.  Such  a  method  is  a  numerical  procedure 
called  the  conjugate  gradient  algorithm  [10].  It  possess  the  same  numerical  stability  as  the 
steepest  decent  algorithm  but  converges  in  a  finite  number  of  steps  at  the  cost  of  some 
additional  computation.  This  method  produces  a  finite  sequence  of  vectors  go,gi,  •  ■  ■  ,<7w, 
where  g n  is  the  solution  we  seek.  The  complete  algorithm  is  given  as  follows.  Let  g0  be 
any  vector  in  RN ,  and  define  do  =  —e0  =  b  —  Qgo,  then  for  A:  =  0, . . . ,  (IV  —  1), 


fffc+i 

dk+ 1 

A 

e* 


gk 

eldk 

~d“Qdk' 

~Cfc+l  +  Pkdlt, 

ef.  ,Q4. 

dlQdt  ’ 

Qs*  -  b. 


(137) 


In  practice,  we  generally  cannot  measure  y(t)  exactly  but  are  given  ym(t)  =  y(t)+n(t) 
where  n(t)  is  a  noise  process.  This,  in  turn,  means  that  we  do  not  know  y  but  are  given 
Vm  =  y  +  n,  where  n  is  the  vector  describing  the  deviation  from  the  true  value  of  y 
because  of  n(t).  Therefore,  application  of  the  conjugate  gradient  algorithm  gives  us  g  in 
a  best  least  square  error  sense.  As  we  will  see,  the  presence  of  noise  in  ym  will  cause 
distortion  and  noise  to  appear  in  the  solution  for  g  and  in  the  resulting  reconstructed 
time  series  derived  from  the  wavelet  expansion. 
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Figures  16  through  18  show  the  input  ®(f),  the  noiseless  output  y(t),  and  impulse 
response  g(t),  the  function  we  seek  to  estimate.  Figures  19  through  21  show  their  re¬ 
spective  wavelet  transforms  that  were  generated  by  using  the  Daubechies  wavelet  of  order 
n  —  5.  Furthermore,  all  wavelet  coefficients  utilized  in  the  deconvolution  algorithm  were 
also  based  on  wavelet  expansions  using  the  Daubechies  wavelet  of  order  n  =  5.  Since  the 
vector  expression  in  Eq.  (132)  is  an  approximation  to  the  matrix  expression  in  Eq.  (130), 
it  is  appropriate  that  the  matrix  D  and  vector  y  (or  ym)  are  composed  of  the  significant 
elements  (those  of  largest  magnitude)  of  the  tensor  V  and  matrix  Y  respectively.  This  is 
done  by  choosing  carefully  the  ranges  of  i,  j,  n,  and  m. 
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Fig  20  —  Square  magnitude  of  the  wavelet  transform  of  the  noiseless 
output  signal  The  transform  kernel  is  a  Daubechies  orthogonal  wavelet 
of  order  n  =  5 


36 


NRL  REPORT  9.116 


Fig.  21  --  Square  magnitude  of  the  wavelet  transform  of  the  true 
impulse  response.  The  transform  kernel  is  a  Daubechies  orthogonal 
wavelet  of  order  n  —  5. 


To  choose  the  ranges  of  i  and  n,  we  consider  the  vector  ym  (or  ym)  whose  elements  are 
the  wavelet  expansion  coefficients  of  the  (measured)  output.  These  coefficients,  in  turn,  are 
equal  to  the  values  of  the  wavelet  transform  of  the  (measured)  output  at  (s,  r)  =  (2J,  2~^m), 
for  j,m  E  Z  x  Z .  Therefore,  we  choose  ji,  jh,  mi,  and  m/,  so  that  y  is  composed  of 
all  coefficients  within  and  on  the  boundary  of  the  region  s  E  [2-2,?1]  and  t  E  [12,26], 
including  the  coefficient  y_ 2,2>  which  is  equal  to  the  wavelet  transform  at  (a,r)  =  (2~2,8). 
As  can  be  seen  from  Fig.  20,  this  region  defines  the  portion  of  the  (.s,t)  plane  where  the 
wavelet  transform  of  the  output  is  significant. 

We  are  now  left  with  the  task  of  choosing  the  ranges  of  i  and  n,  and  this  is  done  by 
examining  the  elements  of  the  matrix  D.  If  we  consider  Eq.  (127),  we  can  assume 
that  the  wavelet  expansion  coefficients  of  any  x  *  <f>(-)  should  at  most  be  significant  over 
the  ranges  of  i  and  n  for  which  the  wavelet  expansion  coefficients  of  the  input  x(t)  are 
significant.  These  coefficients,  in  turn,  are  equal  to  the  values  of  the  wavelet  transform 
of  the  input  at  (s,t)  =  (2*,2~*n),  for  i,n  E  Z  x  Z.  Therefore,  we  choose  ij,  ih ,  nj,  and 
rih  such  that  D  is  composed  of  all  coefficients  within  and  on  the  boundary  of  the  region 
a  d  [ 2 “ 2 , 2 1  ]  and  r  E  [2,16],  including  the  coefficient  d^2,o  that  is  equal  to  the  wavelet 
transform  at  (s,t)  =  (2_2,0).  As  can  be  seen  from  Fig.  19,  this  region  defines  the  portion 
of  the  (s,r)  plane  where  the  wavelet  transform  of  the  input  is  significant. 

The  choice  of  the  ranges  of  i,  j,  n ,  and  m  for  the  example  presented  here  resulted  in 
a  matrix  D  of  size  57x57.  Consequently,  the  conjugate  gradient  algorithm  converged  in 
57  iterations.  The  choice  of  g0  (the  initial  guess  of  g)  is  arbitrary,  hence,  it  was  set  equal 
to  a  vector  of  zeros. 


37 


DAVID  M.  DRUMHELLER 


Figures  22  to  29  show  the  outputs  wdth  various  signal-to-noise  ratios  and  the  resulting 
estimated  impulse  responses.  In  all  cases  the  estimated  impulse  response  was  derived 
from  a  wavelet  reconstruction  using  dilations  2-2,  2_1,  2°,  and  21.  Clearly,  as  the  signal- 
to-noise  ratio  of  the  measured  output  decreases  the  estimated  impulse  response  becomes 
progressively  more  noisy  and  distorted.  This  also  shows  that  high  signal-to-noise  ratios 
are  required  for  accurate  estimation  of  the  impulse  response.  Such  behavior  is  common  to 
all  deconvolution  algorithms;  it  demonstrates  a  typical  trade-off.  Generally,  convolution 
smears  the  impulse  response  to  produce  the  output  signal.  Deconvolution  buys  back  the 
resolution  or  features  of  the  impulse  response  but  at  the  expense  of  producing  an  estimate 
exhibiting  a  signal-to-noise  ratio  that  is  lower  than  the  measured  output  [11]. 
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Fig.  22  —  Noisy  output  signal  with  a  peak-signal- 
to-average-noise  ratio  of  40  dB 
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Fig.  23  —  Estimated  impulse  response  based  on  the  input 
signal  shown  in  Fig.  22  possessing  a  peak-signal-to- 
average-noise  ratio  of  40  dB.  The  wavelet  expansion  for 
the  reconstruction  of  the  impulse  response  used  a 
Daubechies  orthogonal  wavelet  of  order  n  =  5.  and 
dilations  2  " 2 ,  2 _  1 ,  2° .  and  2 1 . 
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Fig.  24  —  Noisy  output  signal  with  a  peak-signal- 
to-average-noise  ratio  of  35  dB 
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Fig.  25  —  Estimated  impulse  response  based  on  the  input 
signal  shown  in  Fig.  24  that  possesses  a  peak-signal -to- 
average-noise  ratio  of  35  dB.  The  wavelet  expansion  for 
the  reconstruction  of  the  impulse  response  used  a 
Daubechies  orthogonal  wavelet  of  order  n  =  5.  and 
di/ations  2  2  \  2°,  and  2l . 


38 


NRL  REPORT  9316 


4.0  1  0  6.0  11.0  16.0  21  0  26  0 

TIME  t 

Fig.  26  —  Noisy  output  signal  with  a  peak-signal- 
to-average-noise  ratio  of  30  dB 


Fig.  28  —  Noisy  output  signal  with  a  peak-signal- 
to-average-noise  ratio  of  20  dB 
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Fig.  27  —  Estimated  impulse  response  based  on  the  input 
signal  shown  in  Fig.  26  that  possesses  a  peak-signal-to- 
average-noise  ratio  of  30  dB.  The  wavelet  expansion  for 
the  reconstruction  of  the  impulse  response  used  a 
Daubechies  orthogonal  wavelet  of  order  n  =  5,  and 
dilations  2-2,  2" 1 ,  2°,  and  21 . 
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Fig.  29  —  Estimated  impulse  response  based  on  the  input 
signal  shown  in  Fig.  28  possessing  a  peal-signal-to- 
average-noise  ratio  of  20  dB.  The  wavelet  expansion  for 
the  reconstruction  of  the  impulse  response  used  a 
Daubechies  orthogonal  wavelet  of  order  n  =  5,  and 
dilations  2~2,  2_1 ,  2°,  and  21 . 


9.  SUMMARY 

This  report  addressed  several  theoretical  and  practical  aspects  of  the  wavelet  trans¬ 
form  and  wavelet  expansion  in  the  context  of  signal  theory  and  signal  processing. 

On  the  most  general  level,  several  theorems  were  proved.  In  particular,  a  ‘decay 
rate  theorem’  was  proved  (Theorem  2)  which  described  how  rapidly  the  wavelet  transform 
decays  as  the  dilation  variable  s  increases.  Moreover,  the  theorem  showed  the  decay  rate 
depends  upon  the  continuity  of  the  transformed  signal.  Such  a  theorem  is  analogous  to 
the  various  decay  rate  theorems  found  in  Fourier  analysis  that  describe  how  rapidly  a 
Fourier  spectrum  decays  as  the  magnitude  of  the  frequency  variable  increases. 
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We  also  presented  the  reformulation  and  extension  of  some  existing  results.  The 
material  in  Section  4  showed  that  linear  system  theory,  i.e.,  input /output  relationships 
for  linear  systems  could  be  reformulated  in  the  context  of  the  wavelet  transform.  It  was 
also  shown  that  stochastic  signal  theory  could  be  applied  to  the  wavelet  transform;  the 
power  spectral  density  and  autocorrelation  function  could  be  used  to  describe  the  expected 
value  of  wavelet  transform  of  a  stochastic  signal.  Some  of  the  material  in  Sections  6  and 
7  extended  some  earlier  work.  Here,  the  continuity,  boundedness,  and  regularity  of  the 
Daubechies  orthonormal  wavelets  were  guaranteed.  Such  results  were  only  sketched  out 
in  the  original  presentation  of  her  work  [2]. 

Two  practical  applications  of  the  wavelet  expansion  were  presented.  The  first  ap¬ 
plication  was  a  filtering  method,  which  may  be  most  useful  when  the  coefficients  of  the 
wavelet  expansion  are  already  available.  Here,  we  showed  that  one  could  passband  filter 
in  the  wavelet  domain.  Because  of  the  sense  of  locality  offered  by  the  wavelet  expansion 
(and  wavelet  transform),  this  filtering  method  may  be  applicable  when  we  require  short 
term,  localized  filtering  of  a  signal.  The  second  application  of  the  wavelet  expansion  was 
the  development  of  a  deconvolution  or  iterative  restoration  algorithm.  The  method  cast 
the  problem  as  a  quadratic  least  squares  problem,  thus  admitting  to  a  solution  by  a  host  of 
well  known  and  established  algorithms.  In  this  case  we  chose  the  conjugate  gradient  algo¬ 
rithm,  because  it  converges  in  a  finite  number  of  iterations.  This  also  allowed  us  to  avoid 
the  problem  of  division  by  zero  that  crops  up  in  the  simple  spectral  division  approach 
to  deconvolution.  The  disadvantage  of  using  the  wavelet  expansion  approach  to  decon¬ 
volution  is  the  need  to  precalculate  the  expansion  coefficients,  and,  to  date,  no  known 
analog  to  the  fast  Fourier  transform  (FFT)  exists  for  the  wavelet  expansion.  Thus,  we 
have  encountered  a  classic  trade  off:  the  development  of  a  robust  deconvolution  algorithm 
at  the  expense  of  additional  numerical  computation. 
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FOURIER  SERIES  COEFFICIENTS  FOR  DAUBECHIES  WAVELETS 


Note  that  the  coefficients  listed  below  define  the  Fourier  series 

H»(f)  =  2~1/2  S'  Ml) 

k= 0 


which  is  consistent  with  Daubechies  definition  of  the  series  Hn(f )  [2].  To  be  consistent 
with  the  form  of  Hn(f)  used  in  Eq.  (82),  one  must  take  hn(k)  =  an(k)/y/ 2. 


a2(0)  =  0.482962913145 
a2(l)  =  0.836516303738 
a2(2)  =  0.224143868042 
a2(3)  =  -.129409522551 


a3(0)  =  0.332670552950 

o3(3)  =  -.135011020010 

a3(l)  =  0.806891509311 

o3(4)  =  -.085441273882 

a3(2)  =  0.459877502118 

a3(5)  =  0.035226291882 

a4(0)  =  0.230377813309 

a4(4) 

=  -.187034811719 

a4(l)  =  0.714846570553 

a4(5) 

=  0.030841381836 

a4(2)  =  0.630880767930 

a4(6) 

=  0.032883011667 

a4(3)  =  -.027983769417 

a4(7) 

=  -.010597401785 

a6(0)  =  0.160102397974 

as(5) 

=  -.032244869585 

a5(l)  =  0.603829269797 

a  5(6) 

=  0.077571493840 

a5(2)  =  0.724308528438 

a«(7) 

=  -.006241490213 

a6(3)  =  0.138428145901 

«i(8) 

=  -.012580751999 

05(4)  =  -.242294887066 

05(9) 

=  0.003335725285 
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ae(0) 

=  0.111540743350 

ae(6) 

=  0.097501605587 

a6(l) 

=  0.494623890398 

a5(7) 

=  0.027522865530 

ae(2) 

=  0.751133908021 

o6(8) 

=  -.031582039318 

a6(  3) 

=  0.315250351709 

ae(9) 

=  0.000553842201 

a6(4) 

=  -.226264693965 

a6(l0) 

=  0.004777257511 

ae(5) 

=  -.129766867567 

a6(ll) 

=  -.001077301085 

<17(0) 

=  0.077852054085 

a7( 7)  =  0.080612609151 

a7(l) 

=  0.396539319482 

ar(8)  =  -.038029936935 

a7(2) 

=  0.729132090846 

a7(9)  =  -.016574541631 

a7(3) 

=  0.469782287405 

a7(10)  =  0.012550998556 

a7(4) 

=  -.143906003929 

a7(ll)  =  0.000429577973 

a7(5) 

=  -.224036184994 

a7(l2)  =  -.001801640704 

a7(6) 

=  0.071309219267 

a7(l3)  =  0.000353713800 

a8(0)  =  0.054415842243 

a*(8)  =  -.017369301002 

a«(l)  =  0.312871590914 

a8(9)  =  -.044088253931 

c8(2)  =  0.675630736297 

a8(10)  =  0.013981027917 

a8(3)  =  0.585354683654 

o8(ll)  =  0.008746094047 

a8( 4)  =  -.015829105256 

a8(12)  =  -.004870352993 

a8(5)  =  -.284015542962 

a8(13)  =  -.000391740373 

a„(6)  =  0.000472484574 

o8(14)  =  0.000675449406 

a8(7)  =  0.128747426620 

a8( 15)  =  -.000117476784 
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a9(0) 

=  0.038077947364 

a9(9) 

=  -.067632829061 

og(l) 

=  0.243834674613 

09(10) 

=  0.000250947115 

*g(2) 

=  0.604823123690 

09(H) 

=  0.022361662124 

a9(3) 

=  0.657288078051 

ag(l2) 

=  -.004723204758 

a9(4) 

=  0.133197385825 

09(13) 

=  -.004281503682 

a9(5) 

=  -.293273783279 

09(14) 

=  0.001847646883 

ag(6) 

=  -.096840783223 

09(15) 

=  0.000230385764 

0.(7) 

=  0.148540749338 

09(16) 

=  -.000251963189 

a9(8) 

=  0.030725681479 

a9(17) 

=  0.000039347320 

<Zjo(0) 

= 

0. 026670057901 

*io(10) 

= 

-.029457536822 

*io(l) 

= 

0. 188176800078 

aio(ll) 

= 

0.033212674059 

*io(2) 

= 

0.527201188932 

*io(12) 

= 

0.003606553567 

aio(3) 

= 

0.688459039454 

*io(13) 

= 

-.010733175483 

*io(4) 

= 

0.281172343661 

Oio(l4) 

= 

0.001395351747 

a10(5) 

= 

-.249846424327 

*io(15) 

— 

0.001992405295 

a10(6) 

= 

-.195946274377 

*io(16) 

= 

-.000685856695 

*io(7) 

= 

0.127369340336 

*io(17) 

= 

-.000116466855 

a10(8) 

0.093057364604 

*io(18) 

= 

0.000093588670 

*io(9) 

— 

-.071394147166 

*io(19) 

= 

-.000013264203 

*11(6) 

= 

0.018692339500 

*n(ll) 

= 

0.031336714900 

*n(l) 

= 

0.144048360129 

*ii(12) 

= 

0.020839548328 

*n(2) 

= 

0.449822419238 

*n(13) 

= 

-.015365977170 

*n(3) 

= 

0.685506451221 

*n(14) 

= 

-.003339972936 

*n(4) 

= 

0.411710892303 

*u(15) 

= 

0.004928945867 

*n(5) 

= 

-.162485521339 

*n(l6) 

= 

-.000308709907 

*n(6) 

= 

-.274320974144 

*n(17) 

= 

-.000893056839 

*u(7) 

= 

0.066025638763 

*n(18) 

= 

0.000249184997 

*n(8) 

= 

0.149791844607 

*n(19) 

= 

0.000054438816 

*n(9) 

= 

-.046504355457 

*n(20) 

= 

-.000034637754 

*n(10) 

- 

-.066445800596 

*n(21) 

= 

0.000004494745 
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<112(0)  = 

0.013114280902 

<112(12) 

= 

0.041627451082 

<112(1)  = 

0.109587064387 

<ii2(13) 

- 

-.012180151045 

<112(2)  = 

0.377449392844 

<ii2(14) 

= 

-.012829445168 

<112(3)  — 

0.657445006413 

<ii2(15) 

= 

0.006713258423 

<112(4)  = 

0.516294170295 

<112(16) 

= 

0.002249393038 

<112(5)  = 

-.044313624533 

ai2(17) 

V 

-.002179176553 

<112(6)  = 

-  .315809615475 

012(18) 

= 

0.000006459278 

<112(7)  = 

-.023471399498 

012(19) 

= 

0.000388621871 

<112(8)  = 

0.182806918672 

012(20) 

= 

-.000088486615 

<112(9)  = 

0.005686977952 

012(21) 

= 

-.000024241195 

<112(10)  = 

-.096186633657 

012(22) 

= 

0.000012775434 

<112(11)  = 

0.010995853244 

012(23) 

= 

-.000001528836 

013(0) 

= 

0.009204916897 

013(13) 

= 

0.002363616024 

013(1) 

= 

0.082889405900 

013(H) 

= 

-.023833745174 

013(2) 

= 

0.312115898739 

013(15) 

= 

0.003917927648 

on(3) 

= 

0.611313131287 

a13(16) 

= 

0.007254616037 

013(4) 

= 

0.589096065406 

013(17) 

= 

-.002760408506 

013(5) 

= 

0.086639694877 

013(18) 

= 

-.001315670455 

013(6) 

= 

-.316237370186 

013(19) 

= 

0.000932006061 

013(7) 

= 

-.126430468961 

013(20) 

= 

0.000049301053 

013(8) 

= 

0.177816118862 

013(21) 

= 

-.000165090932 

013(9) 

= 

0.071915527849 

013(22) 

= 

0.000030664729 

oi3(10) 

= 

-.106342427892 

013(23) 

= 

0.000010440501 

013(H) 

= 

-.026758244166 

013(24) 

= 

-.000004699171 

013(12) 

= 

0.056034390582 

013(25) 

= 

0.000000521846 

46 


NRL  REPORT  9316 


a14(0) 

- 

0.006547491642 

<*i4(14) 

= 

-.029754599557 

a14(l) 

— 

0.063360170581 

<*i4(15) 

= 

-.005754062318 

0.14(2) 

= 

0.259953209778 

<*h(16) 

= 

0.012711190182 

aa4(3) 

= 

0.569486757657 

<*14(17) 

= 

-.000664409841 

ai4(4) 

= 

0.659765991407 

<*14(18) 

= 

-.003831834380 

o14(5) 

= 

0.253248224211 

<*i4(19) 

= 

0.001038385046 

<*14(6) 

= 

-.245883485949 

<*14(20) 

= 

0.000708200880 

«u(7) 

= 

-.207221475070 

<*14(21) 

= 

-.000381870689 

014(8) 

= 

0.141972692112 

<*14(22) 

= 

-.000042656957 

014(9) 

= 

0.144030955893 

<*14(23) 

= 

0.000068164760 

014(10) 

= 

-.083519992219 

<*14(24) 

= 

-.000010124883 

oi4(ll) 

= 

-.071278880702 

<*14(25) 

= 

-.000004370468 

014(12) 

- 

0.054864716315 

<*14(26) 

= 

0.000001706613 

o14(13) 

= 

0.027555092282 

<*i4(27) 

= 

-.000000176357 

